Yang- Yang Anomalies and Coexistence Diameters: Simulation of Asymmetric Fluids 



Young C. Kim 

Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 2074-2 

(Dated: February 2, 2008) 

A general method for estimating the Yang- Yang ratio, IZ^, of a model fluid via Monte Carlo simu- 
lations is presented on the basis of data for a hard-core square-well (HCSW) fluid and the restricted 
primitive model (RPM) electrolyte. The isothermal minima of Ql = {m 2 )\ / (m 4 ) l are evaluated 
at T c in an L x L x L box where m = p — (p) l is the density fluctuation. The "complete" finite- 
size scaling theory for the Q^i n (Tc', L) incorporates pressure mixing in the scaling fields, thereby 
allowing for a Yang- Yang anomaly. It yields a dominant term in the asymmetry, Ql, in — Q^in' 
varying as L~^^ u with an amplitude proportional to the crucial pressure-mixing coefficient, fa. The 
reliably known critical order-parameter distribution for (d = 3) Ising systems then enables one to 
estimate fa, thereby yielding 7£ M , from the Q-minima together with information on the nonuniversal 
amplitudes for the order-parameter and the susceptibility. The detailed analysis needed to estimate 
fa for an HCSW fluid and the RPM is presented. Furthermore, the Q-minima below T c can also 
provide the coexistence- curve diameters, Pdiam(^) = o (P + + P~)> ver y close to T c for both models: 
here p ± (T) are the densities of the coexisting liquid and gas phases. The recently developed recur- 
sive scaling algorithm for Apoo(T) = p + — p~ is adapted to investigate the corresponding universal 
scaling functions. The two extremal forms of these scaling functions are computed with the aid of 
the exactly soluble decorated lattice-gas model. The critical densities for the RPM and HCSW fluid 
found via this route are consistent with previous estimates obtained from the data above T c ; the 
magnitudes of the leading \T — T c \ 213 and \T — r c | 1_c * corrections to Pdiam(T') are estimated. 

PACS numbers: 64.70.Fx, 64.60.Fr, 05.70.Jk 



I. INTRODUCTION 

The presence of a Yang- Yang anomaly means that the 
chemical potential, p a {T), on a phase boundary separat- 
ing liquid and vapor exhibits a divergence of the second 
temperature derivative, d 2 p a (T) / dT 2 , when the critical 
temperature T c of the fluid is approached from below in 
the same fashion as the constant- volume specific heat, 
Cy(T; p c ). This possibility was first proposed by Yang 
and Yang [1] 40 years ago on the basis of a simple ther- 
modynamic relation they derived (referred to later as the 
Yang- Yang equation) which connects the specific heat 
to the pressure and chemical potential derivatives. It 
was only more recently, however, that the Yang- Yang 
anomaly in fluid criticality was seriously investigated 
by Fisher and coworkers [2, 3]. They analyzed care- 
fully experimental data for the two-phase heat-capacity 
of propane (C3H8) and CO2 in the critical region and 
showed that d 2 p<j/dT 2 indeed diverges like the specific 
heat at criticality. Further study of the experimental 
data for propane showed that impurities in the system 
can have significant effects on the heat-capacity data [4] ; 
however, the existence of Yang- Yang anomalies was not 
ruled out. Nonetheless, further careful, experimental in- 
vestigations are desirable [5]. 

Recently computer simulations have become an effi- 
cient tool to study the behavior of fluids and complex 
fluids and, in particular, to enhance our understand- 
ing of their critical behavior. Investigating a Yang- Yang 
anomaly in fluid criticality can thus be aided by simu- 
lations for model fluids, since impurities are absent in 
such models; but such simulations pose a serious chal- 



lenge. Specifically, any sharp divergence of a thermody- 
namic quantity will be rounded owing to finite-size ef- 
fects. These arise from the divergence of the correlation 
length at criticality. In fact, Orkoulas, Fisher and Pana- 
giotopoulos [6] performed grand canonical Monte Carlo 
simulations on a hard-core square-well (HCSW) model 
fluid and concluded that this model probably exhibits a 
negative but small Yang- Yang anomaly quantified by a 
Yang- Yang ratio (see below) IZ^ — — 0.08±0.12; however, 
they could not rule out the possibility that a Yang- Yang 
anomaly was absent. Thus, one might ask: How might 
one measure a Yang- Yang anomaly precisely from simu- 
lations on model fluids? 

To account for a Yang- Yang anomaly, it is necessary 
to mix the pressure deviation, p — p c , into the asymptotic 
scaling fields. More specifically, the full thermodynam- 
ics of a one-component fluid near the bulk critical point 
(p c ,T c ,/i c ) can be described with three relevant scaling 
fields [2, 7], namely, 



p = p- k t- loft + 
t = t-hpb- jip + ■ 
h = fi - kit - j 2 p + 



(1) 
(2) 
(3) 



where the dimensionless deviations of the thermody- 
namic fields from their critical values are 
T - T n . p - p c 



t = 



P ; 



Pck^Tc 



knT r 



(4) 



while ji , J2 , ko, ■ ■ ■ are nonuniversal mixing coefficients 
and p c is the critical density. The scaling hypothesis 
then asserts [2, 7] 

Q\t\ 2 - a W ± (Uh/\t\ A ), 



P 



for t ^ 0, 



(5) 
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where a is the critical exponent for the specific heat while 
A is the gap exponent satisfying the relation A = (J + 7. 
Here j3 and 7 are the critical exponents for the or- 
der parameter and the susceptibility/compressibility, re- 
spectively. In this formulation W±{x) represents two 
branches of a universal scaling function while Q and U 
are nonuniversal amplitudes depending on microscopic 
details of the system under study. Note, however, that 
we have neglected here, for simplicity, both corrections 
to scaling and higher order nonlinear mixing terms [7, 8] . 

The strength of a Yang- Yang anomaly is conve- 
niently measured by TZ tll the limiting ratio of C^iT) = 
—T(d 2 p a /dT 2 ) to the constant-volume specific heat, 
Cv(T, p = p c ) (where p is the number density). It then 
follows that the Yang- Yang ratio IZ^ is related [7, 8] to 
the pressure-mixing coefficient j 2 via 

K, = - ]2 /(l -j 2 ). (6) 

To measure the Yang- Yang ratio, TZ^, quantitatively 
via simulations, the first question is: What thermody- 
namic or other quantities should be studied? The answer 
is not obvious since a direct attempt [6] to study C^T), 
etc. proves not effective. Here we show that the desired 
information can be obtained by carefully investigating 
the finite-size-system parameter Ql defined by [9-11] 

Ql{T; (p) L ) = (m 2 ) 2 L / (m 4 ) L , m = p-(p) L , (7) 

where (-)l is the grand canonical average at fixed T and 
p chosen to yield the desired mean density. In the ther- 
modynamic limit (X — > 00), the parameter Ql then ex- 
hibits surprising, singular behavior on the two sides of 
the coexistence curve [12-14]; namely it vanishes on the 
coexistence curve boundary p = p + (T) and p = p~(T), 
where p^ (T) are the densities of the coexisting liquid and 
gas states in the two-phase region. However, in the one- 
phase region, Q x remains | but drops discontinuously to 
zero at the phase boundary when the two-phase region is 
entered; it then rises continuously to unity as the mean 
density approaches the coexistence-curve diameter [14], 
namely, 

Pdiam(T) = i(p++p-). (8) 

Of course in a finite system, these discontinuities be- 
come rounded. Thus Ql({p)l) exhibits two isother- 
mal minima of magnitudes, say Q m i n (T; L), at densities, 
/o min (T; L), near the coexistence curve boundary [14, 15]. 
(See also Fig. 1 below.) 

For symmetric systems such as Ising models or lat- 
tice gas models, the two minima, <3 min (T; L), have equal 
height while their corresponding densities, /? min (T;L), 
are located symmetrically about the critical density p — 
p c . However, as soon as mixing enters in the scaling fields 
t and h [see (2) and (3)], one finds that Ql((p)l) becomes 
asymmetric. In fact, according to the complete finite- 
size scaling theory (an extension of (5) to finite systems 
[14]), the minima, Q min , and their locations, p min , ex- 
hibit leading asymmetric contributions arising from the 



pressure- mixing coefficient j 2 . More specifically, let us 
define normalized asymmetry factors of the minima via 

A mhl (T;L) = (Q+ n - Q min )/(Q+ n + Q min ), (9) 

B min (T;L) EE Pmin+ + P°»°~_ 2pdiam . (10) 
P min P min 

It is obvious that these quantities vanish identically for 
symmetric systems. For asymmetric situations, however, 
both exhibit leading finite-size correction terms varying 
as L~ fi l v which are proportional to j 2 , while a nonvan- 
ishing (li + ji) combination induces a further correction 
term decaying as L~( A_1 )/ Iy : see below in Sec. III. Here 
v is the correlation-length exponent. Note that for the 
(d = 3)-dimensional Ising universality class to which fluid 
criticality is believed to belong, we have (3/v ~ 0.517 and 
(A - l)/v ~ 0.897: see Table I. It will be shown in Sec. 



TABLE I: Preferred values for the universal Ising d = 3 crit- 
ical exponents adopted here [26]. 
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V 




0.109 


0.326 


1.239 


4.80 


1.565 


0.630 


0.52 



From Ref. [28]. 



Ill that the amplitudes of these leading terms can be 
expressed in terms of the nonuniversal order-parameter 
and susceptibility amplitudes (which, in principle, can be 
measured from simulations) and certain universal con- 
stants. This then raises the hope of measuring j 2 (or 
IZfj) via simulations. 

Indeed, it is shown below that the critical order- 
parameter distribution function for Ising systems and 
the estimated nonuniversal amplitudes, B and C + , for 
the density discontinuity and the susceptibility provide 
precise estimates for the pressure- mixing coefficient, j 2 , 
and, thereby via (6) for the Yang- Yang ratio, TZ^. One 
should notice, however, that measuring j 2 requires pre- 
cise knowledge of T c and p c which must also be found 
by well designed simulations utilizing finite-size scaling 
theory [11, 14]. 

Among the consequences of pressure mixing in the scal- 
ing fields (beyond the divergence of d 2 p a /dT 2 at critical- 
ity) are the appearance of a complex spectrum of singular 
correction exponents in various thermodynamic proper- 
ties [7, 8]. In particular, the pressure-mixing coefficient 
j 2 induces a leading singular term varying as |i| 2 ^ in the 
coexistence-curve diameter, pdiam(T) [7]. This term then 
dominates the previously known |i| 1 ~ Q! term [16]. Hence 
the nonvanishing of j 2 may affect the shape of the diame- 
ter strongly near criticality, and this inevitably hampers 
the precise estimation of the critical density. 

A conventional way of estimating the coexisting liquid 
and gas densities, p^(T), is to observe the grand canoni- 
cal equilibrium distribution function, Vl{p), of the den- 
sity, p, below the critical temperature, T c . In a finite 
geometry of dimensions L d with periodic boundary con- 
ditions the distribution function generally exhibits two 
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peaks near p^(T) when L 3> a where a measures the 
size of particles [9, 15, 17]. For sufficiently large L, these 
peaks can be represented by two Gaussians clearly dis- 
tinguished from each other. Observing these peaks then 
provides direct estimates of p^ (T) when T is sufficiently 
low; for fluids which exhibit no "obvious symmetry" the 
equal- weight prescription has been widely accepted to es- 
timate p 1 (T) from Vl{p) [15, 17]. This approach has 
yielded reasonable estimates of p ± (T) for the HCSW 
fluid [6] and for the RPM electrolyte [18] but only up 
to 1-2% below T c . However, when T approaches T c more 
closely, finite-size effects blur the distinction between the 
coexisting liquid and gas states for most computation- 
ally accessible system sizes thereby seriously hampering 
the reliable estimation of the coexistence curve. Fur- 
thermore, Vl(p) at criticality for finite L still exhibits 
two well-separated maxima; on the other hand, p ± (T) 
should (for /3 > 0) merge precisely at the critical den- 
sity, p c . Hence, unless one can simulate sufficiently large 
system sizes which exceed the correlation length, it is al- 
most impossible to obtain reliable estimates of p ± (T) in 
the critical region via this route; to do so would require 
a rather full understanding of finite-size effects on the 
distribution function. 

To meet this challenge, a scaling algorithm has been 
developed recently for estimating the liquid and gas coex- 
isting densities, p ± (T), of model fluids from grand canon- 
ical simulation data [11, 19]. The algorithm utilizes data 

WminC^ 1 ; L )-> Pmin^; L )) for tne Q-minima and derives 
the density discontinuity, 

Ap oc (T) = (p+ -p~), (11) 

by constructing an appropriate universal finite-size scal- 
ing function. By this route precise estimates for Ap^ (T) 
for the HCSW fluid and the RPM electrolyte were 
obtained up to temperatures only 0.01-0.1% below T c 
[11, 19]. On the other hand, estimating the diameter, 
Pdiam(T), turns out to be more complicated as will be 
demonstrated in the second part of this article. 

In order to calculate the coexistence-curve diameter, 
we will compare the ratio of B m i n to Anin — see (9) 
and (10) — to the average of the Q-minima, Q^ in (T; L). 
As mentioned above, -4 m i n and £> m ; n vanish identically 
in the absence of the mixing coefficients, making the ra- 
tio ill-defined. On the other hand, when the system is 
asymmetric the ratio exhibits rather complex finite-size 
corrections in which both A m in and B m i n contain terms 
varying as L^ 13 ^, i _ ( A_1 )/ I/ , etc. One realizes, there- 
fore, that there is no universal finite-size scaling function 
which yields the diameter in a unique fashion — in con- 
trast to the case of the density discontinuity Ap 00 (T) 
[11, 19]. However, we demonstrate here that there are 
two extremal or limiting cases which yield effective uni- 
versal scaling functions: The first arises when the pres- 
sure mixing is absent or so small as to be negligible: 
see (68) below. The HCSW fluid belongs to this cat- 
egory as demonstrated in Sec. VII A. The second case 
is found when the pressure-mixing term dominates over 



other contributions. An example will be the RPM elec- 
trolyte which then provides us with the other limiting 
scaling function: see (67) and Sec. VII B. Furthermore, 
as for the density discontinuity A Poo (T) [19], these limit- 
ing scaling functions are analytically represented in (104) 
and (107). 

The balance of this article is organized as follows: In 
Sec. II the scaling behavior of the Ql({p)l) is presented 
on the basis of finite-size scaling theories. In Sec. Ill 
the analysis of the Q-minima is presented. Section IV 
is devoted to estimation of the Yang- Yang ratios for the 
HCSW fluid and the RPM by using the critical order- 
parameter distribution of the Ising model and estimated 
nonuniversal amplitudes. Our best estimates for TZ^ are 
given in Eqs. (92) and (95). In Sec. V the scaling al- 
gorithm for estimating the coexistence-curve diameter is 
presented. In Sec. VI an exactly soluble decorated lattice 
gas model is considered to construct one limiting case of 
the universal scaling function for the diameter. In Sec. 
VII we estimate the diameters of the HCSW fluid and 
the RPM electrolyte and compare them with previous re- 
sults. Section VIII summarizes the article and provides 
a brief discussion. 



II. FINITE-SIZE SCALING THEORIES 

In this section we derive the scaling behavior of 
Ql(T, ( P )l) in order to provide the necessary theoreti- 
cal background for estimating the Yang- Yang ratio, IZ^, 
and the coexistence-curve diameter, pdia,m(T). First, 
we consider the large L behavior where one can ob- 
tain Ql(T, (p)l) explicitly on the basis of the double- 
Gaussian approximation. We then study Ql{T, (p)l) in 
the critical region below T c via the complete finite-size 
scaling theory [14]. 



A. Double-Gaussian limit 

For sufficiently large L below T c , we may study 
Ql(T, (p)l) in the two-phase region on the basis of the 
two-Gaussian approximation for the probability distribu- 
tion, Vl{p;T, p), of the density p (at a fixed T and p). 
This may be written 

V L { P] T,p) « C L { X Z 1/2 cxp[-(p-p-fL d /2k B Tx-} 
+ X + 1/2 cM-(p-P + ) 2 L d /2k B T X+ }} 
xexp[p(p-p a )L d /k B T], (12) 

where Cl(p,T) is a normalization constant while the 
X± {T) are the infinite- volume susceptibilities [defined via 
X = (dp/dp)r] at p = p ± (T)±. It is then straightfor- 
ward to calculate the mean density, (p)l, and various 
moments, (m k )L, with m = p — (p)l, for k = 2, 3, • • • , as 
functions of T and p. 
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Following Ref. [14], let us introduce, for simplicity, the 
average and difference susceptibilities 

X(T) = KX++X-) and Xo (T) = |(x+-X-), (13) 

and a convenient parameter (related to the ordering field) 

T(T, h; L) = tanh[% + \xoh)L d ], (14) 

where po(T) = |Apoo(T) [see (11)] and the reduced 
field/chemical potential is defined via 



= [ M - M^)] AbT. 



(15) 



The reduced dimensionless average and difference suscep- 
tibilities are then defined by 



X(L,T) ee x(T)k B T/p 2 L d , 
X (L,T) ee Xo (T)fc B T/^L d . 



(16) 



Note that X and Xo approach zero when L — > oo since 
X and xo ar e finite below T c . 

After integrating (12) multiplied by p, one may write 
the mean density as 



(p) L = Pd iam (T) + xh+ [p (T) + X oh]T. 



(17) 



The expression for Ql can then be written as a ratio of 
two polynomials of fourth order in (p) l as presented in 
Ref. [14]: sec (4. 19)- (4.25). 

In the thermodynamic limit (L — > oo), Ql(T, (p)l) 
becomes zero at the coexistence curve boundary, T = ±1. 
For large but finite L, however, we expect Ql(T, (p)l) 
to exhibit two isothermal minima close to zero at 7± = 
±1 + A7±. To find AT± in terms of X and Xo, we use 
(17) to expand Ql in terms of x oc X and % oc X up 
to linear order to obtain 



[2AT± ± X + X 



4AT±[2 + 3(Xt1 ) + ■■■]" 



(18) 



Solving the equation (dQ l / dT) = for A7± then yields 
the minima at 



AT ± =±\(X ± Xo) + 2 , 



(19) 



where 2 represents terms of order X 1 Xq with i+j > 2. 
On substituting this into (18), one finds 



Q 



± 

min 



X±Xo + 2 . 



(20) 



Note that for the symmetric case (i.e., Xo = 0) these re- 
sults agree with those given in Ref. [14] up to linear order 
in X. For the corresponding densities of these minima, 
we follow [11] and [19] and define the reduced density 
deviation by 



y(T; p) ee 2[p - p Aiam (T)]/ A Poo (T), 



(21) 



where Pdiam(T') is the diameter while Ap 00 (T) is defined 
in (11); here and below we use p to represent the mean 



density, (p)l, unless undesirable ambiguity arises. Note 
that y takes the values ±1 at p = p ± {T). Using (17) and 
(19) yields, after some algebra, 



S/min = ±[l + \XMA/eX)±\XoHA/e 2 X)]- 



(22) 



By taking the mean in (20) and the difference in (22), 
we obtain the scaling relation 



Ay 



mm — 2(^111 ^min) ~ ^ + 2 ^ + 



(23) 



where 



q = Qmin ln[4/eQ min ] with Q min = i(Q+ in + Q min ) 



Similarly, from (9), (10), (20) and (22) one obtains 



(24) 



q = Qminln[4/e 2 Q min ] = q - Q„ 



(25) 



Notice that these relations are universal up to the lead- 
ing order in q (or q) and are asymptotically exact when 
Qmin -> (or L -> oo). 



B. Finite-size scaling closer to criticality 

To understand the behavior of Ql (T, p) near criticality, 
we employ the complete scaling theory (l)-(5) which has 
been extended recently to finite systems of dimensions 
L d with periodic boundary conditions [14]. The finite- 
size scaling ansatz asserts [14, 20] 

Pep w L- d Y(x,z), x = D L tL 1/v , z = U L hL^ /v , 

(26) 

in which Y(x, z) is the basic scaling function while we 
have imposed the hyperscaling relation dv = 2 — a (valid 
for d < 4) and, for simplicity, neglected corrections to 
scaling. Here v is the critical exponent for the correlation 
length, while Dl and Ul are nonuniversal amplitudes of 
dimensions L~ x l v and L~ A//i/ , respectively, which depend 
on the details of the system under consideration: see be- 
low. The scaling function Y(x, z) is universal and thus 
independent of microscopic details of the system; but it 
depends on the geometry and the boundary conditions 
imposed. According to the underlying symmetry, Y(x, z) 
is an even function of z. 

Since the finite-size scaling function, Y(x, z), is ana- 
lytic and even in z, one can expand it as 

Y(x,z) = Y 00 + Y 10 x + Y 2 ox 2 + Y 30 x 3 + --- 

+ z 2 (Y 02 + Y 12 x + Y 22 x 2 + Y 32 x 3 + ■■■) 
+ z i (Y 04 + Y 14 x + Y 24 x 2 + Y 34 x 3 + ■■■) 
+ ■■■. (27) 

One may further normalize Y(x, z) by choosing the 
nonuniversal amplitudes, Dl and Ul, so that one has 
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Yio = Yq2 = 1. The condition Y02 = 1, in particular, will 
be used later for estimating the Yang- Yang ratio, TZ^. 

To recover the bulk limit of the scaling form (5), one 
may set L = l/\Di 1 t\ v and formally let L — > oo. This 
yields (5) with the identifications 

Q=\D L \ 2 - a / Pc and U = U L /\D L \ A , (28) 
W±(z) = Y(±l,z). (29) 

In the two-phase region (t < 0), the scaling function 
W-(z) then has the expansion 

W^{z)^W^ + W^ 1 \z\ + W^ 2 z 2 + ■■■ ■ (30) 

This ensures the density discontinuity at the phase 
boundary. 

The reduced dimensionless number density and suscep- 
tibility are conveniently defined by 

P = P/Pc = {dp/dfi) t , xnn = (d 2 p/dp 2 ) t . (31) 

Similarly one may define the generalized susceptibilities, 
X N k = (d k p/dp k ), k = 3, 4, • • • , which will be used in 
the following section. 



C. Finite-size scaling behavior of Ql 



Finally, the fourth-order susceptibility is 

Xm = e\ P lA^+^l»[{dtY) 

-b ]2 A L-{(d z Y){dtY) + 2(d 2 z Y)(d 3 z Y)} 
-AA l L-\d x dlY) + ---]. (37) 

The desired scaling form for in the grand canonical 
representation (presented in Ref. [11] without any deriva- 
tion) is then readily obtained from (36) and (37). After 
some algebra, we find the central result 

Ql(T, P ;L) = Q Q {x,z)[l+j 2 A j L- K Q j (x,z) 

+A l L- x Q l (x,z) + ■■■}, (38) 

where the universal scaling functions are given explicitly 
by 

Q Q (x,z) = (dlYf/{dtY), (39) 
Q 3 {x,z) = -{d z Y) + lQ{dlY){dlY)/{dtY), (40) 
Ql{x,z) = 4[(d x d?Y)/(dtY) - (d z d x Y)/(d 2 Y)](4l) 

Hence, the symmetry of Y(x, z) implies that Qq(x,z) 
is even in z while Qj(x,z) and Qi(x, z) are odd. The 
pressure-mixing coefficient j 2 thus yields the leading an- 
tisymmetric i-dependent correction term with a decay 
exponent k = (3/v. 



To obtain Ql(T,p) in terms of the scaling variables 
x cx tL 1 / 1 ' and z oc hL A / v , we may notice that Ql can 
be expressed in terms of the generalized susceptibilities, 
XN k , namely, Ql is equivalent to p c V(xnn) 2 /xn 4 - where 
V is the volume of the system. These susceptibilities then 
require the derivatives of the pressure, p, with respect to 
the chemical potential, p, at fixed T. 

In order to compute the susceptibilities in terms of the 
scaling function Y(x, z) in (26), we may, first, obtain the 
reduced density p by differentiating (26) with respect to 
p at fixed T. Using the scaling fields (l)-(3) yields [14] 

p = \ + e x A 3 L-%d z Y) - 32 A,L- K (d z Y) 2 

-AiL~\d x Y) + •••], (32) 

where, setting Iq = 1 [7], the exponents and nonuniversal 
amplitudes are [21] 

ex = 1-J2, k = 0/v, A=(A-l)/i/, (33) 
Aj = Ul/pc, M = (h+h)DL/e 1 U L , (34) 

while, for simplicity, we have adopted the notations 

(d x Y) = {dY/dx), (d 2 x Y) = (d 2 Y/dx 2 ), etc. (35) 

Differentiating (32) with respect to p then yields the 
susceptibility xwjv as 

Xnn = e 2 p c A 2 L^"[(d 2 Y)-3j 2 A 1 L- K (d z Y)(d 2 Y) 
-2A l L~ x {d x d z Y) + -..]. (36) 



D. Finite-size scaling behavior of the reduced 
density deviation, y 

To derive the scaling form for the scale- free density de- 
viation y(T; L), as defined in (21), wc first note that the 
infinite- volume half-discontinuity, po(T) = iApoo(T), 
and the coexistence-curve diameter vary asymptotically 
as 

P0 (T)/p c - B \tf + ---, (42) 
Pdiam(T)/p c = 1 + A 2/3 \t\ 2fJ + A^t] 1 -" + ■ ■ ■ (43) 

where the amplitudes are explicitly given by 

B a = eiQ[/lU-i|T| /3 , A 2(i = -j 2 B 2 / ei , 
Ai_ Q = (2-a)(h+j 1 )QW- \T\ 1 - a 1 (44) 

in which we have introduced the mixing factor [7] 

t = 1 - Mi - (fco + + hh)/ei, (45) 

while Q and U are related to XJ L and D L via the rela- 
tions presented in (28), and W-o and W-\ are expansion 
coefficients of W-(z): see (30). Note that, for simplic- 
ity, we have neglected many higher order terms including 
corrections to scaling [7] . 

On using (32), (42)-(44), one finds, after some algebra, 
that the reduced density deviation (21) has the expansion 

y(T,p;L) - y(x, z)[l + j 2 A j L- K y j (x, z) 

+A l L- x y l {x,z) + -..], (46) 
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where the universal auxiliary scaling functions are 

y(x,z) = (d z Y) /W-i\x\P , (47) 



yj(x,z) = - 
yi(x,z) = - 



{d z Y) + W\\x\W/(d z Y), (48) 
{d x Y) + {2-a)W^\x\ 1 - a 



(d z Y) 



(49) 



Notice that in contrast to the scaling functions, Qq, Qj 
and Qi in (39)-(41), y(x,z), yj(x,z) and yi(x,z) are 
singular at x = 0. Furthermore, y(x,z) diverges as .x|~' 3 
when x — > 0. 

So far we have derived the scaling behavior of 
Ql(T,p;L) and y(T,p;L) in terms of x oc and 
z oc hL A /" near criticality. However, our goal is to ob- 
tain the Q-minima and their locations, Q m [ D (T; L) and 
p^ in (T; L), and to derive scaling relations between them. 
These expressions will then serve to estimate j 2 quanti- 
tatively and, likewise, the diameter near criticality. 



III. ANALYSIS OF THE Q-MINIMA 
A. Q-minima and their locations 

From the observations of simulation data [11, 14, 15] 
and as borne out in the thermodynamic limiting form, 
Ql(T,p;L) is expected to exhibit two isothermal min- 
ima near the coexistence-curve boundary. For symmet- 
ric cases (i.e., l\ = ji = j 2 = 0), it is then obvious that 
Qq(x,z) in (39) must have equi-height minima at some 
value of z, say ±z m (x), for fixed x. Owing to the an- 
alyticity of the finite-size scaling function Y(x, z), the 
function z m (x) is also analytic in x. 

Now we expand the scaling function Y(x, z) about 
these minima at ±z m (x) and, by virtue of the symme- 
try of Y(x, z), obtain 



y(x,z)=^(i)'fl,( I )[ Z T Zm (i)]» 



(50) 



i=0 



where the a* (a;) are universal expansion coefficients of Y. 
Since Y(x, z) is analytic everywhere for finite L, the ai(x) 
must be also analytic and have the universal expansions 



Oi{x) = a i0 + oa a; + a i2 x 2 + ■ 



(51) 



From (38) and (39), in the symmetric case, we then find 
Q^JT; L) « QSun[l + (2a 2 i - a 41 )x + ■ ■ 



■} = Qm(x), 

(52) 

where Q^ in = 0^0/040 is a universal constant. 

In the presence of the mixing coefficients l\, j\ and j 2 , 
the locations of the minima will be shifted by amounts, 
say Az±(x). To find these shifts, we solve the equation 
{dQh/dz) = for Az± perturbatively. After some alge- 
bra, one obtains 



Az±(x) = foAjL-^ix) + AL-^ix^/boix) 



(53) 



where k, A, Aj and Ai were introduced in (33) and (34) 
while the universal scaling functions are 



bo(x) 



, Q4(x) _ a 6 (x) 
a 2 {x) 04(21) 



03(f) 
a 2 (x) 



h(x) = -9a 2 (x) + 10 (a 3 (x)) /a^x), 



b 2 (x) 



a 2 



a' A 
04 



Oi0 3 



,a 3 a 3 



a 4 a 2 



(54) 
(55) 
(56) 



in which a-(x) = dai/dx. 

Substituting these into (38)-(41) via (50) finally yields 

QnuJT;L) = Q m (x)[l±j 2 A jCj (x)L- K ±A l c l (x)L- x +- ••], 

(57) 

where the auxiliary scaling functions, again universal, are 



Cj{x) 

ci(x) 



—ai(x) + Wa 2 (x)as(x)/a4(x), (58) 
4[a' 3 (x)/a 4 (x) - a'^x) / a 2 {x)}. (59) 



Note that the bj(x) in (54)-(56) do not enter here. The 
basic asymmetry factor, Ami n , defined in (9) is thus given 
explicitly by 

A min (T;L) = j 2 A ]Cj (x)L- K + A lCl (x)L- x + ■ ■ ■ . (60) 

As mentioned in the Introduction, we see that the dom- 
inant contribution to -4 m i n arises from the pressure- 
mixing coefficient j 2 with a decay exponent k = [3 jv. 

On the other hand, the locations of the minima, 
y^ in (T;L), can be obtained similarly by using (46)-(49) 
and (53). This yields 



y ± - 



(T;L) = ±y m (x)[l±j 2 A j d j (x)L- K 
±A l di{x)L~ x + -..], 



where the scaling functions are 

y m (x) - -^/ui/s 

dj(x) 



a^xyixfW-u 

a 2 h Wl, 



1 2/3 



-Oi 



oi6 



01 



, , , a' a 2 b 2 W- \x 

di(x) = 1 —-(2 -a) 

Oi oi0 oi 



l-a 



(61) 

(62) 
(63) 

(64) 



while the bj(x) arc defined in (54)- (56). Note that y m (x) 
diverges at criticality as \x\~@ while dj(x) and di(x) ap- 
proach constants singularly as |x| 2/3 and |x| 1_a , respec- 
tively, when x — > 0. 



B. Universal relations for the Q-minima 

The scaling algorithm [11, 19], which was used to esti- 
mate the density discontinuity, Apoo(T), requires a scal- 
ing relation between Ay m j n and Q mnl - From (57) and 
(61), one has 



3>m(a0 



(65) 
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By formally solving for x in terms of Q m i n using (52) and 
substituting into the first member of (65), one obtains a 
function Ay m i n (<3 m in)- Note that its asymptotic behavior 
for Q m in — ► follows from (23). 

The coexistence-curve diameter, pdiam(T), can be esti- 
mated similarly via the scaling algorithm. For this pur- 
pose, we consider 2? m i n , defined in (25), since £> m ; n con- 
tains the diameter only. One should, however, notice 
is defined only for asymmetric cases. Using 
(57) and (61) then yields 

t, _ 32A j d j (x)L- K + A l d l (x)L- x + ■ ■ ■ 

Vmin ~ j 2 A jCj (x)L-« + Aici(x)L~ x + ■■■■ [bb) 

Unlike Ay m i n , 2? m i n has a complicated structure of finite- 
size terms due to the mixing coefficients. One can clearly 
identify two limiting cases. First, when pressure-mixing 
is dominant, one simply has 

Anin - ^4 + 0(L- X+K ) « ej ( X ). (67) 
Cj(x) 

Eliminating x between this relation and (65) and (52) as 
before then yields an asymptotically universal relation 
between Z) m j n and Q m i n . On the other hand, if pressure- 
mixing is absent or can be neglected, the form (66) re- 
duces to 

Vmin = ^ + ---~e l (x). (68) 
ci{x) 

These results will play the central role in the subsequent 
analysis. Note, however, that when Qmin — > one recap- 
tures the universal relation (25) regardless of the partic- 
ular mixing situation. 

IV. DETECTING YANG- YANG ANOMALIES 

In the previous section, the asymmetry factor, 
^4min(T; L), derived from the Q-minima was analyzed 
and it was shown that the pressure-mixing coefficient j 2 
provides a leading term varying as L~PI V which domi- 
nates over other contributions: sec Table I for the Ising 
values for the universal exponents used in the following 
sections. In particular, at criticality [x oc tL x l v = 0) the 
result (60) leads to 

A c min = jiAjCjL-* + A t ciL- x + • • • , (69) 

where the universal constants are 

Cj = -aio + 10a 2 o»3o/«4o, (70) 
ci = -4an/a 2 o + 4a 3 i/a 40 . (71) 

in which are the coefficients of x^ in the expansion 
of Y(x, z) at the minima, z — z m : see (51) above. On 
the other hand, Aj is a nonuniversal constant introduced 
in (34), while Ai oc (l\ + ji)/(l — jz)- Thus to measure, 
in particular, the pressure- mixing coefficient, ji-, one re- 
quires not only the simulation data for A^ in but also 
information concerning the universal constant Cj and the 
nonuniversal amplitude Aj . 



A. Estimation of the universal constant Cj 

To determine the constant Cj as given by (70), one must 
determine the expansion coefficients, a^o, of the scaling 
function Y(x,z) at x = about its minima: see (50). 
Invoking universality, we thus consider a simple cubic 
ferromagnetic Ising model, to which universality class 
normal fluids are believed to belong. In this case the 
"pressure-field" , p c p, in (26) corresponds to the reduced 
free energy density for a system of volume L d while the 
ordering field h corresponds to the reduced dimension- 
less magnetic field. If we now define a scaled fluctuating 
magnetization density by w = A^-mL 13 ^ (where to corre- 
sponds to the fluctuating reduced magnetization density) 
with the identification Al = l/UL,we find from (26) that 

(w k ) = {{A L mL p/v ) k ) = (d k Y/dz k ). (72) 

Thus the expansion coefficients ajo can be obtained from 
the j-th moment of the scaled magnetization density at 
criticality. 

Now in order to estimate the (w k ) numerically, we will 
utilize the most reliable and precise estimate available for 
the universal probability distribution function, P£(m), of 
the magnetization density, to, at criticality: this has been 
obtained numerically from careful simulations [22] . This 
function may be written in scaling form as 

P£(to) = Cexp{-/(to) +wz}, w = A L mL^" ', 

(73) 

where the nonuniversal amplitude Al is chosen so that 
the distribution has unit variance in w, while C is a nor- 
malization constant; but one should notice from (27) and 
(72) that this is equivalent to the normalization con- 
dition Y 2 = 1 for Y(x,z): see Sec. II B. Recall that 
z = UjJiL^I 1 ' with Ul = 1/Al. Here f(w) corresponds 
to a universal canonical scaling function. Tsypin and 
Blote [22] present this function approximately in the form 

f(w) - [(w/w ) 2 - 1] 2 [a(w/w ) 2 + c] . (74) 

where Wo, a and c are constants estimated via the simu- 
lations as 

w = 1.134(4), a = 0.158(2), c= 0.776(2), (75) 

where the uncertainties in parentheses refer to the last 
decimal place quoted [22] . Furthermore, the condition of 
unit variance yields Al ~ 0.802. 

Using these forms, one may calculate the derivatives of 
Y"(0, z), i.e., (d k Y/dz k ), numerically by computing the k- 
th moments of w. It is then straightforward to calculate 
Ql{T c , (p)l) as a function of z: see Fig. 1. It exhibits 
two symmetric minima, at z = ±z^ where 

Q c nlin = 0.1178(15) and z^ = 2.24(3). (76) 

Note, however, that Q^m ^ s about 7% higher than the 
value (~ 0.1107) found for the HCSW fluid via grand 
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FIG. 1: Plot of the universal critical-point scaling function 
Ql(T c ;( P }l) w Qq(0,z) vs z = U L hL A/v for the (d = 3) 
Ising model obtained numerically from (73) [22]. The dashed 
lines locate the minima at ±Zm ~ ±2.2395. 

canonical simulations [11, 23]. The universal constant, 
Q c , at criticality (i.e., z — 0) agrees with the estimated 
value for the Ising universality class, Q c ~ 0.6236 [10, 24]. 
Calculating (w k ) at the minimum, z = z^, then yields 
the desired expansion coefficients as 

a 10 ~ 1.2213, a 20 ^0.1751, o 30 ~ -0.1271, 

a 40 ~ 0.2603, a 60 ^ 0.9878. (77) 

Notice that in (52) one obtains Q^ in = 020/040 — 0.1178 
which is fully consistent with the estimate (76). Finally, 
the universal constant cj in (70) is found to be 

cj ~ -2.077. (78) 

In reality the distribution of the magnetization, P£(m), 
should behave like exp(—cs\m\ s+1 ) for large to with S = 
A/f3 ~ 4.8 for (d = 3) Ising universality class: see Table 
I. Hence, the approximate form (74) docs not describe 
the behavior of the distribution properly for large to, even 
though § + 1 ~ 5.8 is close to the exponent 6 embodied in 
this form. One can improve the approximation of Tsypin 
and Blote in order to capture the correct behavior of large 
to. However, we find that appropriate modifications leave 
the central estimate (78) for cj unchanged. 

B. Nonuniversal amplitude Ul 

The nonuniversal amplitude U l determining the value 
of Aj via (34) depends on the microscopic details of the 
system. Unlike the universal constant Cj obtained in the 
previous section, one must estimate Ul separately for 
each model. Here we present a method for evaluating 
Ul from knowledge of the leading nonuniversal critical 
amplitudes for bulk thermodynamic quantities. 



In bulk near-critical systems, the susceptibility above 
T c and the order parameter below T c vary as [25, 26] 

M(T) = {dp/ dp) « B\tf, (79) 
Xnn(T) = (d 2 p/dp 2 )^C + r\ (80) 

with p = p/k^T and p, = fi/k^T, where the nonuniversal 
amplitudes can be obtained from the scaling formulation 
(5) as [7] 

B ee Pc B = (l-j 2 ) Pc QUW^\Tf, (81) 
C+ = (1- J2 ) 2 PcQU 2 W +2 \t\-<, (82) 

while t is defined in (45). Note that these amplitudes 
have dimensions L~ d . Recall that W-\ and W +2 arc 
expansion coefficients of the scaling functions W±(z): see 
(29) and (30). Thus they are universal. Likewise note 
that Q and U are related to the finite-size amplitudes, 
D L and U L , via (28). 

Now we may combine (81) and (82) in order to solve 
for Ul in terms of B and C + . After some algebra, one 
finds 

U L = UD£ = K [B~>(C+f] 1/(2 - a) /\l- 32 \, (83) 
where if is a universal constant given by 

/ „ x -1/(2-q) 

K=(w\Wl 2 ) , (84) 

which can be obtained numerically from the well-studied 
Ising model. As mentioned above, we first notice that the 
condition of unit variance and the numerical simulations 
[22] yield U L = 1/A L ~ 1/0.802 for the (d = 3) nearest- 
neighbor simple cubic (sc) Ising model [22]: see (73). The 
amplitudes B and C + for this model have been studied 
and are given by [25, 27] 

5 = 1.66(3), C+ = 1.095(10). (85) 

From (83) and (85) together with the exponent values in 
Table I, one then obtains 

K = 0.881(12). (86) 

Now, one clearly sees from (83) that to estimate Ul 
for model fluids, information concerning the amplitudes 
B and C + is crucial. 



C. Estimation of the Yang- Yang ratios 

In this section we finally obtain estimates for the 
pressure-mixing coefficients j 2 , and thereby the Yang- 
Yang ratios, TZ^, for the HCSW fluid and the RPM by 
utilizing the results obtained. In Fig. 2 the asymmetry 
factors, *4 m i n (L; T c ), at the critical temperature for the 
two models are presented as functions of _L _ ( A_1 )/ Iy and 
L~PI V , respectively. If a Yang- Yang anomaly is present, 
the data should decay asymptotically as L~ /3/V when 
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FIG. 2: Plots of the asymmetry factors, _4Jj lin (Z/), as obtained 
from simulations for (a) the HCSW fluid (vs L _ ( A_1 )/ I/ ) and 
(b) for the RPM (vs L^ 13 ^) with the exponent values listed 
in Table I. The solid curves are fits of the data to (69). 



L —> oo. The plot for the RPM clearly suggests that this 
highly asymmetric model has a nonvanishing Yang- Yang 
anomaly. On the other hand, the HCSW fluid would 
seem to have a quite small, if any, Yang- Yang anomaly. 
The solid lines are fits of the data to the formula (69) 
neglecting the higher order terms. [See (91) and (94) be- 
low.] Note that the fit for the HCSW fluid exhibits a 
small negative leading amplitude. Based on these data 
— the origin of which we first review briefly — and the 
resulting fits we now describe the procedures used to es- 
timate the Yang- Yang ratios, H^, quantitatively. 



1. Details of simulations 

The asymmetry factors presented in Fig. 2 and all the 
simulation data in the subsequent figures have been ob- 
tained via grand canonical Monte Carlo (GCMC) simu- 
lations in a cubic box of volume V — L 3 . In GCMC one 
performs a simulation at a thermodynamic state point 
(SP) characterized by a given value of the temperature, 
T, and the chemical potential, p. In order to capitalize on 
the widely used multi-histogram reweighting technique 
[30], one measures the joint histogram of the fluctuat- 
ing energy U and the number of particles N for different 



SPs [6, 24]. This approach enables one to extract the 
maximum amount of information from the simulations. 

To obtain sufficiently accurate simulation data for the 
HCSW fluid, each simulation was performed for a total 
run length in the range of 20-80xl0 6 MC steps, depend- 
ing on the particular system size under investigation [6]. 
A total of 30-100 SPs (which broadly cover the critical 
region and the coexistence curve) have been used in the 
computations. Statistical uncertainties for the density 
p = N/V and for the leading moments for each histogram 
are found to be less than |%. On the other hand, for the 
RPM the total run length needed for each simulation is 
considerably larger owing to the much lower critical tem- 
perature, and to the great number of SPs necessary to 
ensure accurate data. For example, a total of 167 SPs 
were used for L = 12a (where a is the particle diame- 
ter), in which a typical SP has a (2-10) xlO 4 independent 
samples amounting to ^10 10 MC steps [24]. For smaller 
systems, the number of independent samples employed is 
generally larger by factors of 10-50. For each SP, the sta- 
tistical uncertainties of the raw data are less than 1-2%. 

We know of no definitive or systematic study con- 
cerning the propagation of errors in the multi-histogram 
reweighting process. However, experience shows that 
when the SPs are closely spaced in regions of rapid 
changes in the computed quantities that encompass de- 
sired values of T and p, the systematic uncertainties in 
the leading moments do not exceed and may well be 
appropriately less than the statistical uncertainties that 
characterize each SP. Hence, we believe that the errors 
associated with the data presented in Fig. 2 and the sub- 
sequent figures are no larger than the symbol sizes (or, 
in plots like Fig. 3, no more than a couple of line thick- 
nesses) . Furthermore, the confirmation of a smooth and 
systematic variation of specific plots with increasing val- 
ues of L, as seen in Fig. 3, serves as an important cross- 
check. When apparent erratic behavior is seen, larger 
runs and/or further SPs have been employed as a check. 



2. Hard-core square-well fluid 

The HCSW fluid investigated here consists of hard- 
spheres of diameter a with an attractive square-well of 
depth e in the range a < r < 1.5a, where r is the interpar- 
ticle distance. Previously Orkoulas et al. [6] tried to esti- 
mate Tift via Monte Carlo (MC) simulations by comput- 
ing directly the second derivative of the chemical poten- 
tial along the phase boundary; they concluded that there 
was a negative but small Yang- Yang ratio, 7£ M ~ —0.08. 
Owing to the finite-size effects, however, they could not 
rule out a vanishing 71^ . Here we use the asymmetry fac- 
tor, ^4 m in, at criticality as illustrated in Fig. 2 to calculate 
Hp more precisely. In doing so, as mentioned above, one 
needs to compute Ul which requires knowledge of the or- 
der parameter (or density-half-discontinuity), po(T), and 
the susceptibility, Xnn{T), above T c . 
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Near criticality the appropriate expansions are 

Po (T) = B\tf[l + b e \t\ e + ■■■], (87) 
Xnn(T) = C+t-^l + cgt 6 + ■■■}, (88) 

where 9 is the leading correction-to-scaling exponent. For 
(d = 3) Ising criticality one has 9 ~ 0.52 [28], while 
bg and eg arc nonuniversal amplitudes. Recent precise 
simulation of the coexisting densities for this HCSW fluid 
via the Q-minima scaling algorithm (see below) yields 
a reliable estimate for B: see Fig. 1 in Rcf. [11]. On 
the other hand, to determine C + , we have computed the 
finite-size susceptibility xnn(T; L) by simulations on the 
critical isochore p = p c . In the bulk limit XNN{T)\ty 
should then approach C + when t — ► 0: see Fig. 3(a). For 



0.10 




FIG. 3: Plot for estimating the susceptibility amplitude C 
from data for \nn on the critical isochore [6, 24]. Note that 
t — (T — T c )/T c while the values ip = 1 (solid lines) and 
tp = 9 — 0.52 (dotted lines) have been chosen for the auxil- 
iary extrapolation exponent (together with, purely for clarity 
of presentation, A\ = 1 and A$ = 0.4, 0.8). For (a) the 
HCSW fluid and (b) the RPM electrolyte, the estimated crit- 
ical densities are p* = p c a A = 0.3206 and 0.079, respectively, 
while L* = L/a = 9, 10.5, 12, 13.5 for the HCSW fluid, and 
7, 8, 9, 10, 12 for the RPM. The dashed lines represent approx- 
imate upper and lower extremal estimates for C + . 

finite systems, however, this product will approach zero 
when t — > 0, as seen in Fig. 3, since xnn(T; L) is always 
bounded for L < oo. Nevertheless, we can estimate C + 



reasonably well by extrapolating the data of [11] to t = 0. 
We conclude 

Ba 3 = 0.612 ± 0.005, C + a 3 = 0.076 ± 0.005, (HCSW). 

(89) 

Furthermore, the Ising fit to the po(T) for the HCSW 
fluid via (87) with a further correction term varying as 
|i| 2/3 [7] yields the leading correction-to-scaling amplitude 
as bg ~ -1.04. Using (83) and (86) now yields the HCSW 
nonuniversal amplitude as 

U L 1 1 -h\a^ jv = 0-410 ±0.012. (90) 

The fit of the Anm data to (69), illustrated in Fig. 2, 
yields 

j 2 A j c j ~ -0.117, Aici ~ 0.570, (HCSW), (91) 

where the two amplitudes have opposite signs. Finally by 
using (34), (78), (90) and (91) and the estimate p c a 3 ~ 
0.3067 [6], we obtain 

n = 0.040 ± 0.003, 7^ = -0.042 ± 0.003 (HCSW). 

(92) 

This is, in fact, consistent with the previous, less precise 
estimate — 0.08 ± 0.12 [6]. It should be noted that the 
agreement of the two, quite different procedures for esti- 
mating TZ,j. serves as an encouraging check on the overall 
validity of the scaling analysis. Note that the uncertainty 
quoted in this estimate reflects the uncertainties only in 
B, C + and p c , but not in the fitting coefficients in (91): 
because of the higher order terms it is hard to provide a 
realistic estimate of these uncertainties. It is reasonable 
to believe, however, that they are likely to be no more 
than 20-30%. 

The amplitude Ai is proportional to l\ + ji- see (34). 
One may then hope to estimate l\ +ji, if not separately, 
from the fit (91) via the same route. To do so, however, 
requires further information concerning the nonuniversal 
amplitude Dl and the universal constant q: see (34) and 
(71). The former can be obtained via the normalization 
condition Fio = 1 in the expansion (27), as for Ul\ but 
this requires knowledge of the temperature dependence 
of the finite-size-scaling function Y(x, z). This is also the 
case for estimating the universal constant q which con- 
tains an and 031. Hence, further investigation is needed 
to obtain estimates for these subdominant mixing coeffi- 
cients. 



3. Restricted primitive model 

The restricted primitive model (RPM) electrolyte con- 
sists of an equal number of positive and negative ions, of 
charges ±go and hard-core diameter a interacting with 
each other via the Coulomb potential, ip{r) = ±q^/Dr, 
where D is the dielectric constant of medium. For this 
model the Yang- Yang ratio TZ^ has not previously been 
investigated seriously. The markedly asymmetric nature 
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of the coexistence curve hints that IZ^ might be large 
compared to the HCSW fluid. Indeed, studies of the 
generalized /c-susceptibility-loci [14] have suggested that 
1Z^ might be positive and nonnegligible. 

To estimate 7Z^ for the RPM (at a £ = 5 fine-lattice 
discretization level [24]) we follow the procedure used for 
the HCSW fluid. Thus Fig. 3(b) presents the effective 
susceptibility amplitude XNNt 1 , along the critical iso- 
chore p = p c as obtained from grand canonical MC sim- 
ulations [24]. An Ising fit to the precise data for the 
coexistence curve [11] can also be well achieved. Thence 
we find 

B = 0.284 ± 0.01, C+ = 0.087 ± 0.005, (RPM) (93) 

with b e ~ -1.17 in (87). 

As seen in Fig. 2, grand canonical MC simulations for 
the RPM yield strong asymmetry values. The fit of the 
A^ in {T c ;L) data to (69) provides the amplitudes 

fcAjCj ~ 1.644, Aici ~ 0.395, (RPM). (94) 

It is remarkable that j 2 Aj has the opposite sign to that 
for the HCSW fluid while the amplitude combination is 
also ~ 14 times larger. Together with the critical density 
estimate, p c a 3 ~ 0.079 [24], we thus find 

j 2 = -0.35 ± 0.07, Uf, = 0.26 ± 0.04 (RPM). (95) 

In contrast to the HCSW fluid, IZ^ for the RPM is posi- 
tive and large which seems to reflect the strongly asym- 
metric nature of this model electrolyte. As for the HCSW 
data the uncertainties in (95) do not reflect those in (94); 
however, as the combination jiAj is rather large and the 
fit in Fig. 2 rather good further uncertainties of only, say 
5 to 10%, seem likely from this source. 

One should notice that the RPM studied here is, as 
mentioned above, a discretized version with the dis- 
cretization parameter ( — 5 [24]. It has been shown, 
however, that the universal critical behavior is indepen- 
dent of £ (> 4) [29]. Furthermore, the nonuniversal crit- 
ical parameters, T c and p c , for £ > 5 are close to the 
continuum limits; the T c (( — 5) and p c (( = 5) arc only 
3% and 5% higher than the continuum values, respec- 
tively. Thus we believe that TZ^iC = 5) is likely to be 
quite close to TZ^C, = oo). 

V. SCALING ALGORITHM FOR THE 
COEXISTENCE -CURVE DIAMETER 

Estimating the coexistence-curve diameter and, in par- 
ticular, identifying its singular behavior near criticality 
has been a major challenge for both experiment and sim- 
ulation. As mentioned in the Introduction, the presence 
of a nonvanishing pressure-mixing coefficient 32 yields a 
|i| 2/3 term in the diameter, that dominates the previously 
anticipated [16] |i| 1_a contribution. Here we present a 
scaling algorithm designed to enhance the estimation of 



the diameter near criticality, and, thereby, to improve 
estimates of the critical density, p c . 

The recently developed algorithm for estimating the 
density discontinuity, Ap 00 (T), utilizes the scaling rela- 
tions between the Q-minima and their normalized loca- 
tions, y^ in : see (57) and (61). The asymptotically exact 
expression for large L at fixed T < T c given in (23) pro- 
vides the limiting guide to construct a universal scaling 
function (65) by finding optimal values for Ap^T). The 
step-by-step procedure for implementing the algorithm is 
presented in Refs. [11] and [19]. Here we adapt the al- 
gorithm to derive the coexistence-curve diameter. For 
this purpose we consider the relation between T> min [con- 
taining /9diam(T) as a variable: see (8)-(10) and (25)] and 
Qmin. This exhibits a rather complicated structure of 
finite-size corrections owing to the various mixing coeffi- 
cients: see (66); however, the relation achieves a univer- 
sal form (25) for large enough L. Thus we may hope to 
construct a scaling function from this limiting behavior 
via a scaling algorithm by optimally choosing values for 
/Odiam(T). Although the steps for the Ap OC) (T) algorithm 
are presented in detail in [19], we recapitulate the main 
points in the present setting for the sake of completeness. 

There are three main steps: (i) Collect data 
sets of the minima of Ql and their locations, 
{Qmin( T ; L 0:Pmin( T ;^)}, for a range of box sizes 
{Li}f =1 at fixed values of T < T c . For this purpose, 
multi-histogram rewcighting technique [30] should be em- 
ployed to generate the data at any desired temperature 
for a given system size. In practice, n — 3 distinct box 
sizes with L3 > 1.3Li may well suffice; to avoid the ef- 
fects of the corrections to scaling one finds that one needs 
at least L\/a > 8 (where a measures the size of particles). 
However, this guide will surely depend on the system un- 
der consideration, (ii) Choose a value T = Tq sufficiently 
low that the two-peaked, double-Gaussian structure of 
Vl{To;p) is well realized. In this case Qmin is close to 
zero. (For the HCSW fluid we obtained Q min < 0.03.) 
Now note that in this region the universal scaling rela- 
tion between V min and Q m i n in (25) is exact up to linear 
order. Hence, at this temperature we select a diameter 
estimate, say px , independent of the L t , which leads to 
the best fit of 

^min = {Pm + Pm ~ 2 PT }/ (Pm ~ Pm)/Anin 

VS. q^=q(T ;Li) (96) 

to the relation (25). The value pr can then be identified 
as an estimate for pdiam(To). (iii) Increase T by a small 
amount AT to T\ — T + AT and compute X> m i n (7i; Li) 
and qf 1 = q(Ti; Li). In doing so, it is necessary to choose 
AT small enough that the new set {q[ overlaps the 

previous one {<Zo^}™ =1 . Note that q increases with T, ap- 
proaching the critical value q c at T c . When the new data 
set is obtained, we choose, as before, a new value, pr 1 , 
so that the new data set collapses optimally onto the 
previous one; this procedure then extends the previously 
computed scaling function to larger values of q. Again, 
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the new value pr 1 can be regarded as an estimate for 
/9diam(Xi)- These steps are repeated iteratively by in- 
creasing the temperature Tj to Tj + \ = Tj + ATj. This 
extends the scaling function further and generates succes- 
sive estimates for pdiamiTj) for j = 2, 3, • • • . When criti- 
cality is approached, smaller increments ATj are needed 
and high quality data prove essential. For graphical il- 
lustrations of these procedures, see Fig. 2 of Ref. [19]. 

One should note that owing to the competitive na- 
ture of the singular terms in (66), with k ~ 0.517 and 
A ~ 0.897, there is no leading universal function that 
will, in general, yield the diameter. This is in contrast to 
the case of the density discontinuity where there is a uni- 
versal scaling function which provides an easier and faster 
algorithm to estimate Apoo(T): see [19]. It is therefore 
necessary to approach the problem case by case. In the 
following sections, we study the two extremal limits and 
compute two distinct universal scaling functions: see (67) 
and (68). First, the scaling relation (68) derived in the 
absence of ji can be obtained from an exactly soluble 
decorated lattice gas model which is known to exhibit a 
nonvanishing l\ but has no pressure mixing (i.e., j 2 = 0). 
One may suspect that the HCSW fluid, which has been 
seen to exhibit a rather small pressure-mixing coefficient 
j 2 may be well approximated by this route. The other 
extremal case, in (67), in which pressure mixing appears 
to strongly dominate, can be obtained by analysis of the 
RPM electrolyte. 



VI. DECORATED LATTICE GAS MODEL 

Consider a decorated lattice gas model which can be 
solved exactly in terms of the solution of the associated 
Ising model [31, 32]. The model exhibits a entropy- like 
singularity, in the diameter but has no pressure 

mixing and so serves as a guide to the extremal case (68) 
for the diameter. 

The decorated lattice gas we study here consists of pri- 
mary cells centered on sites of a basic simple cubic (sc) 
lattice and secondary, decorating cells centered on the 
bonds between the nearest neighbor pairs of sc lattice 
sites. All cells have equal volume, say vo, and do not 
overlap. Each cell can be empty or occupied by at most 
one particle. Particles in nearest neighbor primary cells 
then interact with energy — e while particles occupying 
neighboring primary and secondary cells interact with 
energy — Ae. For simplicity we will consider only A = 1 
which will suffice for our present purposes. Of course, 
one may consider a simpler version by only allowing the 
interaction between particles in nearest primary and sec- 
ondary cells. However, this will not change the significant 
results. 

The grand canonical partition function of the model 
can be expressed in terms of that of the ordinary lattice 
gas model [31] . Let N be the number of primary cells and 
q the coordination number (with q = 6 for the sc lattice) 
so that the total number of cells is N to t — {q + 2)N/2. 



The dimensionless activity of molecules in the decorated 
lattice gas is z = ^oA^ 3 exp(p / IcbT) , where At is the de 
Broglie wavelength and p is the chemical potential [33], 
and we write K = e/ksT for the coupling constant of the 
decorated lattice gas model. Then the grand canonical 
partition function can be written as 



E(z,K) = (l + Z y N / 2 E(z,K), 



(97) 



where S is the corresponding partition function of the 
ordinary lattice g function of the transformed ac- 

tivity and coupling constant 



l + ze XK ^ q 



l + z 



K = K + \n 



ze 



2XK\1 



{l + ze XK Y 



(98) 
(99) 



The critical point values, K c and z c , follow by substitut- 
ing K c = 4:Kl s and z c = cxp(— 2gif < ! s ) where is the 
critical coupling of associated sc Ising model. 

The coexistence curve of the decorated gas can be simi- 
larly obtained from that for the Ising model. The number 
density in the decorated gas is given by 

p = lim N t -l(dlnZ/d\nz) 
qz 



(q + 2)(l + z) 



+ 



q + 2 



p(z,K) 



<91nz 
d\nz 



w(z, K) 



dK 
<91nz 



,(100) 



where p and —w are the number density and the energy 
per site of the ordinary lattice gas, which can be written 
on the phase boundary as 



p± = |(1 ± M ) 



w± 



l(y + u±6M ) 



(101) 



in which Mq and u are the spontaneous magnetization 
and the energy density of the Ising model below T c ; here 
± represent the liquid and gas sides of the phase bound- 
ary, respectively. Note that the phase boundary in the 
activity-coupling plane is given by z(K) = cxp(—qK/2). 
Liu and Fisher [25] give approximate forms for the spon- 
taneous magnetization and the specific heat below T c as 

M (T) ~ B\tf{\ -a^l 1 / 2 ), 
C{T) ~ ^|tr Q (l- acff |i| 1 / 2 ) + 6 cff , (102) 

which are valid through the whole temperature range. 
It is worth mentioning that one can calculate l\ exactly, 
while the pressure-mixing coefficients, ji and ji-, are iden- 
tically zero; but computing l\ is not relevant here. 

Now using (100) and (101), and the precise numer- 
ical results for the spontaneous magnetization, Mq(T), 
and specific heat for the Ising model obtained via the 
fitting formulae (102) together with the amplitude val- 
ues in [25], we can obtain the coexistence curve, p ± (T), 
for the decorated lattice gas model via (100). Thus 
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FIG. 4: Scaling plot of (Ay„ 



-1//3 



vs. q 



ln(4/eQ 



(with (3 = 0.326) for the coexistence curve of the decorated 
lattice gas model. The various symbols depict results gener- 
ated from simulations at increasing temperatures with system 
sizes, L* = 8, 9, 10 and 11. For comparison, the universal scal- 
ing curve obtained previously, and represented analytically in 
Eq. (12) of [19], is presented as a dashed curve. 




FIG. 5: Plot of the (subdominant) scaling function ei(x) vs. 

Q — Qmin 

ln(4/e Qmin) for the diameter of the decorated lat- 
tice gas model with A = 1. The dashed line represents the 
two-Gaussian limit (25) which is asymptotically exact to lin- 
ear order in q, while the solid curve is a fit to (104). 



we can construct the scaling function for Ay m ; n pre- 
sented in Fig. 4 as symbols. One should note, however, 
that, although the coexistence curve for the decorated 
gas is exactly known, building the scaling function for 
Ay m i n requires information concerning the Q-minima; 
these are finite-size quantities, which are not known ex- 
actly. Thus, to obtain the Q-minima and their locations, 
{Qmin( T '' L )iPnvn( T '> L )}i for the decorated gas model, 
we have performed grand canonical MC simulations on 
lattices of dimensions L x L x L with periodic bound- 
ary conditions. The resulting scaling curve in Fig. 4 is 
indistinguishable from that of the HCSW fluid [19] - 
see the dashed curve — built via the recursive scaling al- 
gorithm, thus confirming our iterative scaling approach 
for constructing Ap QO (T). We then build a scaling plot 
for e/(x) — see (68) — by using the exact diameter of 
the decorated lattice gas. This curve, presented in Fig. 5, 
should represent the case when j 2 = and l\ ^ 0, namely 
Anin ~ Q/di. We expect that an effective scaling curve 
for the HCSW fluid should be close to this one. 



VII. NUMERICAL ESTIMATES FOR THE 
DIAMETERS 



Following the procedures explained we now report ex- 
plicit results for the HCSW fluid and the RPM. 



A. Hard-core square-well fluid 

To derive the diameter for the HCSW fluid, we start 
from the temperature Tq = k^T/e = 1.10 at which, as 
for the density discontinuity [19], the double-peak Gaus- 
sian for the density distribution is quite accurate. At this 
temperature, we determine that Pdiam(Io) which leads to 
the best fit of £> m ; n (To;L) vs. q(Tq; L) to the asymptot- 
ically exact two-Gaussian limit given in (25): see the 
dashed straight line in Fig. 6(a). Since £> m i n depends 
on asymmetry, high quality data are required as stressed 
before. As seen in Fig. 6(a), the data exhibit relatively 
noisy behavior at low temperatures (possibly owing to 
the lack of enough histograms covering the whole range 
of the two-phase region at this temperature). The scatter 
may be compared to the analysis of the density discon- 
tinuity (see Refs. [11, 19]) where the cancelation of the 
uncertainties on the two sides of the coexistence region, 
when the average is taken, evidently results in smoother 
behavior. Despite the scatter in Fig. 6(a) we can estimate 
the diameter at T to within about ±0.2% relative to p c . 
Following the procedure in Sec. V, we then build up a 
scaling function numerically by increasing the tempera- 
ture, Tj, and finding optimal values for p<nam(Tj). The 
function constructed this way for the HCSW fluid is pre- 
sented in Fig. 6(a); it should approximate the extremal 
scaling limit — see e/(x) in (68) — in which pressure 
mixing is considered negligible. 

The scaling analysis of the Q-minima indicates that 
Anin^) ~ ei(x) = di{x)/ci(x) in (68) converges to a 
universal constant, say Ci, when q approaches the critical 
value q c = q c - Q c miD ~ 0.175 as \q - q c \^ a . From (59), 
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FIG. 6: Plot of scaling functions for the coexistence-curve di- 
ameters vs. q = q — Qmin for (a) the HCSW fluid and (b) the 
RPM. The symbols represent data at different temperatures 
while the dashed lines represent the exact, two-Gaussian lim- 
iting behavior (25). For comparison, the plots display, as solid 
and dotted curves, fits to the extremal scaling curves, ei(x) 
and ej(x), using (104) and (107). 



(64), (68) and (71), this constant C/ is given by 

-l 

C l =c l /d l (0) = c l 



«01 , 020^2 



O-10 Oio^OO 



(103) 



where boo and b 2 o can also be expressed in terms of the 
via (54) and (56). Computing Ci, however, requires val- 
ues for ctoi, an and 031, which can be obtained from the 
temperature dependence of the finite-size scaling function 
Y{x,z); but, at this stage these values are not known. 
Nevertheless, we have estimated C\ via extrapolation 
from the HCSW data and obtained Ci ~ 0.418. Us- 
ing this value, we can then fit the scaling curve to the 
approximant 



ei(x) ~ Ci 



1 



(l- q ')^—)(l + Siq i + 32q ii) 



1 + hq' + t 2 q' [ 



, (104) 



where q' = q/q c while we set t\ = si — 1 + a + q c /2Ci in 
order to ensure the small-g limiting behavior (25). The 
values following from this fit are shown in both parts of 
Fig. 6 as solid curves: the fitted coefficients are si ~ 



4.53, s 2 ^ -6.61, s 3 ~ 1.20, h ~ 3.85, t 2 ^ -9.08, and 
i 3 ~ 4.24. Note the singular behavior at T c when q' — > 1 
implied by the (1 — g') 1_a factor in (104) which yields a 
vertical tangent in the plot. 

Figure 7(a) presents the diameter, Pdiam(r), for the 
HCSW fluid obtained via the scaling algorithm: see solid 
circles. For comparison, previous estimates obtained di- 
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FIG. 7: Plots of the coexistence curve diameters, PdiamC^*) = 
Pdi am a 3 obtained via the scaling algorithm: solid circles, (a) 
The HCSW fluid with T* = k B T/e and (b) the RPM (at a 
fine discretization level ( — 5 [24]) with T* = k^TDa/qi. 
The crosses are previous estimates [6, 18] obtained from data 
for the two-peak structure of the density distribution via an 
equal-weight prescription. The open circles in the insets are 
the estimates of the critical point, (T*,p*). 

rectly from the simulated probability distribution func- 
tion of the density are presented (as crosses) . In the tem- 
perature range where the equal- weight prescription works 
well, the estimates from the two methods agree within 
the uncertainties. On the other hand, the equal-weight 
prescription cannot provide reliable values for the diam- 
eter in the near-critical region. The current approach, 
however, gives Pdia,m(T) far closer to the T c . A fit to the 
asymptotic expansion [7] 

Pdiam(T) ~ Pc [1 + A 2p \t\ 2f} + A^ a \t\ X - a + Alt] , 

(105) 

yields the critical density as p* = p c a 3 = 0.3072 ± 0.0005 
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which agrees with the previous central estimate [6] within 
the uncertainties , while the amplitudes are A 2 p — 0.37, 
Ai- a ~ -2.14, Ax ~ -2.52. Note that the magnitude 
of A 2 p is almost an order smaller than that of Ai- a . 
The sharp curvature for the diameter very close to T c 
reflects the |t| 1_Q singularity. The magnitude seems sur- 
prisingly large but appears as an inescapable conclusion 
of our analysis: its validity might benefit from further 
investigations and, indeed, experiments. 

As discussed above, we anticipate that when T> m [ n (q) 
is constructed from data for the HCSW fluid in moder- 
ately small systems, it will match the scaling form (68); 
of course, when j 2 is not identically zero, it should even- 
tually reveal the limiting form (67) when L — > oo. In fact, 
this is confirmed by comparing Figs. 5 and 6(a), where 
the scaling curve for the HCSW fluid (solid curve) is seen 
to be in good agreement with that of the decorated lat- 
tice gas model. On the other hand, for the RPM a strong 
j_,-&l v term arises from j 2 ; in this case the behavior of 
the X\mn scaling function should approach ej(x) identi- 
fied in the asymptotic form (67). To check this, we move 
on to the calculation of the diameter for the RPM. 



B. Restricted primitive model 

No changes in the previous procedures are needed to 
derive the diameter for the RPM. Using the data set 
of the Q-minima employed in [11, 19], we construct the 
scaling function for T> mln starting from the low temper- 
ature T* = k B T/{q 2 /Da) = 0.0426 ~ 0.84T*. The scal- 
ing function so constructed is presented in Fig. 6(b): it 
clearly differs from that for the HCSW fluid. As ar- 
gued, this is simply due to the relatively large pressure- 
mixing coefficient of the RPM. The scaling analysis yields 
Anin ~ e j( x ) m (67) which approaches a universal con- 
stant, say Cj, as \q — q c \ 2 ^ when T — > T c — . Just as for 
Ci, the constant Ca can be written as 



Cj=c j /d j (0) = c j 
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(106) 



where 6 o = b (0) and b w = fei(0). Using (77) and (78) 
via (54) and (55), we have Cj ~ 0.87. One may then fit 
the scaling curve to the further approximant, as for the 
HCSW fluid, 



e,{x) ~ Cj 



1 - 



(W) 2/3 (1 + -<W + -S2q' 2 ) 



l + tiq 1 + t 2 q' 2 



(107) 



with h = S! - 2/3 + q c /2Cj. We obtain si ~ -1.871, 
s 2 ~ 1.050, s 3 - -0.169, ii ~ -2.422, t 2 ~ 1.952 and 
t$ ~ —0.529. The fit is presented in Fig. 6 as dotted 
curves. As one might expect from the exponent value 
2/3, this plot rises more sharply to the limiting constant 
than does ei{x) with exponent 1 — a. 

In Fig. 7(b) we display the present estimates (solid 
circles) for pdiam(T) for the RPM (at a fine discretization 
level ( — 5 [24]) along with the previous values (crosses) 



[18] estimated from the density distribution via the equal- 
weight prescription. The agreement is good within the 
uncertainties; but, as anticipated, the current approach 
yields reliable estimates much closer to the critical point. 
Extrapolation provides the critical density estimate p* c = 
p c a 3 ~ 0.078(3) in excellent agreement with the previous 
estimate p* ~ 0.0790(25) [24]. Furthermore, a fit to (105) 
yields A 2[j ~ -2.03, A^ a ~ 28.3 and A x ~ -23.5. It is 
interesting that although 2(3 ~ 0.65 is less than 1 — a ~ 
0.89 the numerical behavior of the RPM diameter as T 
approaches T c seems smoother than for the HCSW fluid: 
note, however, that the estimates for the RPM approach 
only to \t\ <~ 10~ 3 whereas those for the HCSW fluid 
show sharp behavior almost a decade closer to T c . 



VIII. SUMMARY 

In summary we have provided a general method for de- 
termining the strength of the Yang- Yang anomaly from 
simulations of model fluids. Specifically, we have stud- 
ied the isothermal minima of the fourth-order fluctua- 
tion parameter, Ql(T;p), in detail on the basis of the 
two-Gaussian approximation, that is exact well below 
T c , and of the complete finite-size scaling theory near 
criticality [14]. It was shown that the asymmetry fac- 
tor, A m - m oc (Q^in _ Qmin)' exhibits a leading term de- 
caying as L~P/ U and of magnitude set by the pressure- 
mixing coefficient, j 2 , followed by a L~( A ~ 1 )/ 1 ' term aris- 
ing from the combination of the mixing coefficients l\ and 
j\: see (60). We then showed that precise finite-size data 
for .4 m i n at T c provide a quantitative route to estimat- 
ing the pressure- mixing coefficient, j 2 , and thereby the 
Yang- Yang ratio, TZ^. By using universal information 
for the critical order-parameter distribution of (d = 3) 
Ising systems and the specific critical amplitudes of the 
order-parameter and the susceptibility above T c for the 
model fluids under study, one can estimate j 2 rather pre- 
cisely. This method was applied to the HCSW fluid with 
range-to-core ratio 1.5 and the RPM electrolyte, leading 
to Uf, ~ -0.042 and +0.26, respectively: see Sec. IV C. 
The approach can be applied readily to any model fluid 
system: it will be a challenge to understand which fea- 
tures of a system govern the sign and magnitude of . 

We have also presented in detail a recursive scaling 
algorithm using the Q-minima which enables one to esti- 
mate precisely the liquid-gas coexisting densities, p^(T), 
very close to T c . Corresponding universal scaling func- 
tions which, in principle, can be derived from the prob- 
ability distribution function, Vl(T; p), were investigated 
numerically via the algorithm by using grand canonical 
simulation data for the HCSW fluid and for the RPM 
[6, 24]. The two leading extremal universal scaling func- 
tions for the diameter were calculated and represented 
analytically in (104) and (107). The algorithm yields 
precise results for pdiam{T) in a range of temperature a 
decade or two closer to T c than was previously feasible: 
see Fig. 7. The new estimates for the critical tempera- 
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tures and densities for both models agree well with the 
best previous estimates extrapolated from the data above 
T c . Furthermore, the behavior of p<uam(T) close to T c for 
the two models compares favorably with experimental 
data for SF 6 and liquid metals, respectively [34] 

The new method is applicable to any model for which 
the order-parameter distribution can be reliably estab- 
lished at temperatures well below the critical tempera- 
ture. To obtain successful estimates, however, one needs 
high quality data for the Q-minima and their locations. 
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